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Abstract 

Several  algorithms  for  introducing  artificial  dissipation  into  a  central  difference  approx¬ 
imation  to  the  Euler  and  Navier  Stolces  equations  are  considered.  The  focus  of  the  paper  is 
on  the  convective  upwind  and  split  pressure  (CUiSP)  scheme,  which  is  designed  to  support 
single  interior  point  discrete  shock  waves.  This  scheme  is  analyzed  and  compared  in  detail 
with  scalar  and  matrix  dissipation  (MATD)  schemes.  Resolution  capabUity  is  determined 
by  solving  subsonic,  transonic,  and  hypersonic  flow  problems.  A  finite-volume  discretiza¬ 
tion  and  a  multistage  time- stepping  scheme  with  multigrid  are  used  to  compute  solutions 
to  the  flow  equations.  Numerical  results  are  also  compared  with  either  theoretical  solutions 
or  experimental  data.  For  transonic  airfoil  flows  the  best  accuracy  on  coarse  meshes  for 
aerodynamic  coefficients  is  obtained  with  a  simple  MATD  scheme. 
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Introduction 


Accuracy  must  be  a  primary  consideration  in  the  construction  of  any  numerical  scheme.  One 
would  like  to  devise  a  scheme  with  the  minimum  amount  of  artificial  dissipation  required 
for  stability,  as  well  as  convergence  in  the  case  of  a  stationary  solution.  For  fluid  dynamic 
computations  the  numerical  scheme  should  be  designed  to  have  hi^  axxuracy  in  smooth  regions 
of  the  flow  field  and  well  resolved  shock  waves  and  contact  discontinuities.  According  to 
Harten  [3]  such  discrete  formulations,  where  the  accuracy  away  from  discontinuities  is  at  least 
second  order,  are  called  high  resolution  schemes.  The  design  of  these  schemes  for  systems  of 
conservation  laws  is  generally  based  on  theory  developed  for  a  scalar  conservation  law.  As  a 
consequence  one  cannot  ensure  that  the  properties  of  the  scheme  for  the  scalar  equation  are 
valid  for  the  system.  In  addition,  schemes  that  permit  high  definition  of  shock  waves  without 
oscillations  are  first  order  in  the  neighborhood  of  shocks.  Concern  naturally  arises  regarding 
contamination  of  the  solution,  especially  in  the  case  of  viscous  flows.  For  these  reasons  the 
properties  and  resolution  capability  of  this  class  of  schemes  must  be  verified  through  numerical 
applications  for  a  wide  range  of  flow  conditions. 

High  resolution  schemes  of  particular  interest  for  solving  the  compressible  fluid  equations 
are  those  that  allow  shock  capturing  with  a  single  interior  point.  In  Ref.  [6]  Jameson  presents 
the  convective  upwind  and  split  pressure  ( CUSP)  scheme.  For  this  scheme  the  artificial  diffusive 
flux  vector  associated  with  a  given  coordinate  direction  is  expressed  in  terms  of  changes  in  the 
state  and  flux  vectors.  A  somewhat  limited  number  of  inviscid  and  viscous  computations  have 
been  performed  to  evaluate  these  schemes  (see  Refs.  [6,  7]  and  [22,  23,  9]), 

In  the  current  work  we  investigate  and  analyze  the  HCUSP  version  [6]  of  the  CUSP  scheme. 
The  HCUSP  scheme  allows  a  solution  with  constant  total  enthalpy  for  steady  flow.  We  discuss 
the  shock-capturing  behavior  for  various  choices  of  the  dissipation  coefficients.  We  introduce 
a  simple  modification  of  the  limiter  function,  which  is  generaUy  used  with  the  scheme,  to 
control  background  dissipation,  and  thus  global  accuracy.  The  HCUSP  scheme  includes  a 
contribution  that  is  scaled  according  to  the  local  velocity.  If  the  velocity  approaches  zero, 
as  it  does  for  viscous  flow  near  a  solid  boundary,  and  there  is  a  high  aspect  ratio  mesh, 
the  dissipation  in  the  streamwise  direction  (i.e.,  direction  of  long  side  of  mesh  cell)  may  not 
be  adequate  for  convergence.  A  change  in  the  velocity  scaling  factor  based  on  aspect  ratio  is 
presented.  The  resolution  capabUify  of  the  HCUSP  scheme  is  evaluated  for  subsonic,  transonic, 
and  hypersonic  flow  problems.  A  detailed  comparison  of  the  scheme  with  a  scalar  and  matrix 
dissipation  schemes  is  performed.  The  scalar  scheme  is  based  on  the  dissipation  model  of 
Jameson,  Schmidt,  and  Turkel  [4]  and  is  often  used  with  central  differencing. 

Dissipation 

A  finite-volume  approach  is  applied  to  discretize  the  fluid  dynamic  equations  of  motion.  The 
computational  domain  is  divided  into  quadrilateral  cells,  fixed  in  time,  and  for  each  cell  the 
governing  equations  can  be  nondimensionalized  and  written  in  integral  form  as  follows: 

^  [[  wdxdy+  [  {fdy  -  gdx)  =  f  {f^dy  -  g^dx),  (1) 

otJJn  Jm  Re  Jdn 

where  H  is  a  generic  cell  with  boundary  911.  In  the  scaling  factor  for  the  viscous  terms  on  the 
right  hand  side  of  Eq.  (1),  the  quantities  7,  M,  and  Re  are  the  ratio  of  specific  heats,  Mach 
number,  and  Reynolds  number,  respectively,  with  M  and  Re  defined  in  terms  of  nominal  con¬ 
ditions.  Taking  Wj^k  as  the  ceU-averaged  solution  vector,  Eq.  (1)  can  be  written  in  semidiscrete 
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form  as 


where  ilj^k  is  the  area  of  the  cell,  and  £  is  defined  by 

C  =  Cc  +  Cd  + 

with  the  subscripts  C,  D,  and  AD  referring  to  convection,  diffusion,  and  artificial  dissipation. 
In  order  to  simplify  the  description  of  the  dissipation  models,  we  consider  the  one-dimensional 
Euler  equations  of  gas  dynamics. 

Scalar  Dissipation  Model 

The  scalar  dissipation  is  based  on  the  model  introduced  by  Jameson,  Schmidt,  and  Turkel 
[4].  This  model  defines  a  switching  function  based  on  a  blending  of  the  second  and  fourth 
differences.  The  term  associated  with  the  operator  Cad  is  expressed  as 

CAUWj  =  -{D'^  —  D^)wj  =  (ij+i/2  —  dj  1/2- 

Then 

DS-  =  V[(A,.+V2£fi/2)A]^,-,  (2) 

DS-  =  v[(A,.+i/2  45i/2)AVA]u;,-,  (3) 

where  the  index  j  refers  to  a  cell  center,  and  the  operators  A  and  V  are  forward  and  backward 
difference  operators*  The  variable  scaling  factor  Aj4-i/2  is  defined  as 

^j+i/1  =  2[('^f)j  +  (■^f)j+i]> 

where  is  the  largest  eigenvalue  in  absolute  value  (i.e.,  spectral  radius)  of  the  flux  Jaco¬ 
bian  matrix  associated  with  the  the  Euler  equations.  In  two  dimensions,  with  denoting 
arbitrary  curvilinear  coordinates,  the  scaling  factor  is  usually  defined  as 

{h)iM  =  =  1  +  '^j,k 

with  r  =  The  exponent  C  is  generally  taken  to  be  between  ^  and  |.  The  coefiicients 

and  use  the  pressure  as  a  sensor  for  sharp  gradients,  and  they  are  defined  as 

4+1/2  =  max(i/j_i,  Vj,  i/j+i,  Vj+2), 

Pj  i+^Pj+Pj+i 

4?i/2  =  [°’  “  4+1/2)]  ’ 

where  typical  values  for  the  constants  and  are  in  the  ranges  3  to  ^  and  ^  to 
respectively.  We  shall  refer  to  Eq.  (4)  as  the  JST  switch.  The  switching  function  v  can  be 
interpreted  as  a  limiter,  in  the  sense  that  it  maximizes  the  second-difference  contribution  at 
extrema  and  switches  off  the  fourth-difference  term.  Moreover,  at  shock  waves  the  dissipation 
is  first  order,  and  a  first-order  upwind  scheme  is  produced  for  a  scalar  equation.  In  smooth 
regions  of  the  flow  field  the  dissipation  is  third  order. 
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Thus,  we  have  two  different  dissipation  mechanisms  at  work.  The  switch  determines  which 
one  is  active  in  any  given  r^ion.  For  smooth  flows,  i/  is  small  and  the  dissipation  terms  consist 
of  a  linear  fourth  difference  that  damps  the  high  frequencies  which  the  central  difference  scheme 
does  not  damp.  This  is  necessary  for  achieving  a  steady  state.  In  the  neighborhood  of  large 
gradients  in  the  pressure,  i/  becomes  large  and  switches  on  the  second-difference  viscosity 
while  simultaneously  reducing  the  fourth-difference  dissipation.  This  is  needed  to  introduce 
an  entropy  condition  so  that  the  correct  shock  relations  are  satisfied  and  to  reduce  overshoots 
near  discontinuities. 

Matrix  Dissipation  (MATD) 

Sharp  resolution  of  shock  waves  without  oscillations  can  be  achieved  by  closely  imitating  an 
upwind  scheme  in  the  neighborhood  of  a  shock  wave.  A  key  feature  of  upwind  schemes  is  a 
matrix  evaluation  of  the  numerical  dissipation.  With  this  evaluation  the  dissipative  terms  of 
PHfh  discrete  equation  are  scaled  by  the  appropriate  eigenvalues  of  the  flux  Jacobian  matrix 
rather  than  by  the  spectral  radius,  as  in  the  JST  scheme.  A  matrix  dissipation  model  can 
easily  be  constructed  by  starting  with  the  JST  formulation. 

One  can  show  [20]  that  the  necessary  modification  of  the  JST  scheme  to  produce  a  matrix 
dissipation  model  is  the  substitution  of  |A|  for  the  eigenvalue  scaling  factor  A  in  Eqs.  (2)  and 
(3).  Since  the  Euler  equations  are  a  strongly  hyperbolic  system,  the  coefficient  matrix  can  be 
diagonalized.  Assume  QAQ  ^  =  A  (diagonal  matrix).  Then  |A|  is  defined  as  |A|  =  Q  ^|A|Q 
and  [A]  =  dm^dAil,  IA2I,  [Aal),  where  Aj  are  the  forward  acoustic,  backward  acoustic,  and 
convective  eigenvalues.  Efficient  ways  of  computing  |A|  times  a  vector  are  presented  in  Ref.  [20]. 

In  practice  one  cannot  choose  Ai,A2,A3  as  the  eigenvalues.  Near  stagnation  points  A3 
approaches  zero  while  near  sonic  lines  Ai  or  A2  approaches  zero.  A  zero  artificial  viscosity 
would  create  numerical  difficulties.  Hence,  we  limit  these  values  as 

[Ajl  =max(|Aj|,y„p(A)),  j  =  l,2 

IA3I  =max(|A3|,14p(A)), 

where  p{A)  is  the  spectral  radius  of  A,  and  the  linear  eigenvalue  A3  can  be  limited  differently 
than  the  nonlinear  eigenvalues.  The  parameters  I4,  have  been  determined  by  numerical 

experimentation.  Typical  values  are  Vn  =  0.25  and  Vi  =  0.025. 


TVD  and  LED  Properties 

In  one  dimension  consider  the  approximation 

=  Wj  -  ^(/j+1  -  /i-l)  +  ^(Qj+l/2^’"i+l/2  -  Qj  1/2), 

where  =  y^j+\  ~  Suppose  that  the  Jacobian  matrix  A  is  at  least  locally 

constant.  Then  the  scheme  is  TVD  if 

Qj+1/2  ^  |^J+l/2|  •  (^) 

This  result  follows  directly  from  the  fact  that  the  system  can  be  diagonalized,  and  we  can 
ensure  that  each  characteristic  variable  satisfies  the  TVD  criterion  (i.e.,  the  total  variation  of 
the  variable  is  nonincreasing). 
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A  replacement  for  the  switch  of  Eq.  (4)  that  is  TVD  for  a  scalar  equation  is  introduced  in 
Ref.  [20].  In  one  dimension  this  switch  is  given  by 

^  ^  bj'+i  ~  +Pj-il  ^0^ 

^  \pj+i-pj\  +  \pj  -Pj  il  +  £ 

and  choose  =  ^.  In  practice  we  usually  use  a  weaker  form  than  Eq.  (6),  for  example, 

\Pj+i-‘^P}+Pj  i| 

^  (1  —  Loypj’VD  “I" 

where 


Ptvd  =  \Pi+\  -Pjl  +  \Pi  -Pj  il> 

P  =  Pj+i  +  ^Pj  +  Pj-ii 

and  0  <  a;  <  1.  The  TVD  switch  of  Eq.  (6)  is  recovered  when  w  <  1.  Typically  oj  ~  1/2. 

In  Refe.  [6]  and  [8]  Jameson  develops  a  theory  for  scalar  nonoscillatory  schemes  based  on 
the  local  extremum  diminishing  (LED)  principle  that  maxima  do  not  increase  and  minima  do 
not  decrease.  The  LED  principle  applies  to  multidimensional  problems,  and  it  is  equivalent  to 
the  TVD  principle  in  one  dimension.  If  a  scheme  is  LED,  then  the  scheme  is  positive. 


SLIP  Scheme 


Until  now  we  have  considered  a  combination  of  a  low-order  and  high-order  artificial  viscosity 
based  on  a  scalar  switch.  This  switch  has  the  disadvantage  that  it  is  based  on  only  one  quantity, 
the  pressure.  Moreover,  it  forces  all  variables  to  be  treated  equal,  even  though  some  experience 
sharp  changes  throu^  the  discontinuity  while  others  are  continuous  across  the  shock.  One 
can  instead  limit  independently  each  dependent  variable  in  each  coordinate  direction. 

Jameson  [5]  introduced  a  new  class  of  limiters  and  implemented  such  a  limiting  process 
within  a  framework  that  he  called  the  symmetric  limited  positive  (SLIP)  scheme.  For  the  SLIP 
formulation  Jameson  constructed  a  family  of  limiter  functions  based  on  the  function 


R(u,  v)  =  1  — 


u  —  v  ^ 
\u\  +  |i;|  +  e 


where  q\a  a.  positive  number  and  e  has  the  dimensions  of  u.  Note  that  R(u,  u)  «  0  whenever  u 
and  V  have  the  opposite  sign.  Let  w  be  an  element  of  the  solution  vector  for  the  flow  equations. 
According  to  our  previous  theory  [20]  R(Awj^-i/2i  where  Awj^s/2  =  '^j+2  —  u'j+i, 

would  be  replaced  by  i'j+i/2)  where  Vj+1/2  is  the  maximum  of  i/j  over  the  nearest  neighbors 
and  V  is  given  by  Ekj.  (6).  Define  the  limiter  function  L(u,v)  by 

L(u,  v)  —  R(u,  v)  ■  (7) 


At  the  mesh  cell  interface  j + 1/2,  we  define  the  left  and  right  states  for  each  dependent  variable 
as 

Wi  =  Wj  +  ^L(AWj+ii2,AWj_ii2), 

1 


(8) 


and  so 


Wr  —  'Wl  =  ~  ^(^'“^j+3/2)  ^'^’j  1/2)- 

For  the  artificial  viscosity  all  differences  will  be  based  on  wr  —  wi.  In  the  nei^borhood  of 
shock  waves  R{u,v)  and  hence  L{u,v)  are  close  to  zero.  Moreover,  wr  —  wr  =  ^Wj^i/21 
resulting  in  a  first-order  scheme  for  the  artificial  viscosity.  Fbr  smooth  flow  R[u,  v)  ~  1,  and 
L{u,v)  =  {u  +  v)/2.  Hence,  in  a  smooth  region 


WR  -WL  =  ^Wj+i/2 
~  AWj+i/2 


L{Awj+z/2,^Wj  1/2) 
^'^’j+Z/2  +  ^Wj_i/2 

2 


(9) 


Thus,  the  SLIP  scheme  has  similar  properties  to  the  JST  scheme.  One  can  obtain  the  rela¬ 
tionship  between  the  SLIP  and  JST  schemes  by  defining  the  diffusive  flux  dj+1/2  as 


<^^•+1/2  =  Olj+l/zbl’R-'^L) 


(10) 


with  ctj^i/2  =  K^^^-^j+i/2-  The  quantity  is  a  parameter,  and  A  is  the  spectral  radius  of  the 
associated  flux  Jacobian  matrix* 

One  difference  between  the  JST  and  SLIP  schemes  involves  the  parameters  and 
for  the  second  and  fourth  differences,  respectively.  Both  and  are  free  parameters  in 
the  JST  scheme.  As  seen  from  Eqs.  (10)  and  (9)  these  parameters  are  automatically  chosen  as 
and  with  the  SLIP  scheme.  The  coefficient  of  the  second  difference  is  chosen  as  5 
so  that  the  scheme  is  fully  upwind  for  supersonic  flows.  However,  the  purpose  of  the  fourth- 
difference  viscosity  is  only  to  accelerate  the  convergence  to  a  steady  state  by  eliminating  the 
odd-even  point  decoupling.  Hence,  we  wish  to  be  as  small  as  possible  for  accuracy  while  stUl 
achieving  a  good  convergence  rate.  It  does  not  seem  reasonable  to  connect  the  two  components 
of  the  artificial  viscosity. 

One  can  generalize  the  SLIP  scheme  by  reintroducing  a  free  parameter  that  essentially 
governs  the  level  of  the  third-order  viscosity  in  smooth  regions.  The  resulting  scheme  has  the 
disadvantage  that  a  free  parameter  must  be  chosen;  however,  it  has  the  advantages  of  greater 
flexibility  and  increased  accuracy.  We  repleice  Eq.  (7)  by 


L(u,v,w)  =  R{u,v)]  • 


(1  - 


(11) 


where  left  and  right  state  values  are  determined  by 

WL  =  Wj  +  ^L(A'ii;j+3/2,  AiOj+i/2,  Aii;_,_i/2), 
WR=Wj+i  -  iL(A'u;j+3/2,Au;j.(,i/2,Aii;j  1/2)- 


When  =  1,  the  L  of  Eq.  (11)  reduces  to  the  original  L  of  Eq.  (7).  Now,  at  shock  waves 
R{u,w)  is  small  and  we  have  wr  -  wl  =  Awj+1/2.  For  smooth  regions  of  the  flow  field  we 
have  WR  —  wl  =  —2K^‘^^A!^Wj^i/2- 
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CUSP/AUSM  Scheme 

So  far  we  have  described  the  use  of  an  artificiaJ  viscosity  based  on  either  a  scalar  or  matrix  co¬ 
efficient.  Liou  and  coworkers  designed  a  scheme  called  Advection  Upstream  Splitting  Method 
{AUSM)  [10, 11, 12].  This  method  was  later  refined  for  large-scale  3-D  viscous  computations 
[15].  AUSM  is  based  on  a  splitting  of  the  flux  function  into  convective  and  pressure  contribu¬ 
tions.  In  some  sense,  the  pressure  terms  contribute  to  the  acoustic  waves  while  the  velocity 
terms  contribute  to  convective  waves.  Hence,  it  is  reasonable  that  these  flux  terms  be  treated 
differently.  Liou  thus  considers  decompositions  of  the  flux  vector  that  are  not  based  on  a  char¬ 
acteristic  decomposition  but  on  Mach  number  scaled  contributions  of  the  left  and  the  right 
states  to  the  interface  flux.  This  decompositton  has  the  disadvantage  that  it  is  more  difficult 
to  develop  for  other  sets  of  equations  compared  with  a  characteristic  decomposition.  A  similar 
type  scheme  called  the  convective  upwind  split  pressure  (CUSP)  scheme  was  later  introduced 
by  Jameson  [5]  and  subsequently  modified  by  Tatsumi,  Martinelli,  and  Jameson  [7,  8,  22,  23]. 
The  CUSP  scheme  has  several  advantages.  First,  one  can  consider  the  scheme  as  another  type 
of  artificial  viscosity,  since  it  is  defined  as  a  sum  of  the  central  flux  average  plus  a  dissipative 
flux.  Hence,  it  can  be  readily  used  with  a  variety  of  time-stepping  schemes  (e.g.,  multistage, 
LU,  implicit,  etc.).  Second,  the  CUSP  formulation  can  be  used  with  multistage  schemes  which 
do  not  evaluate  the  artificial  dissipation  fluxes  at  every  stage,  in  order  to  reduce  computa¬ 
tional  work.  Another  advantage  of  the  CUSP  scheme  is  that  it  can  be  easily  combined  with 
preconditioning,  since  preconditioning  is  based  on  the  inviscid  flux  form  and  not  the  artificial 
dissipation.  Hence,  we  shall  only  describe  the  CUSP  version  of  this  type  of  scheme. 


Definition  of  CUSP  Scheme 

Previously,  we  introduced  the  scalar  and  matrix-valued  viscosities  by  considering  dj^i/2  of  the 
form  ^ 

dj+i/2  =  2*3j+i/20^'’j+i  ~ 

The  factor  1  is  introduced  so  that  we  get  first-order  upwinding  when  Qj+1/2  =  !■  We  note 
that  for  the  scheme  to  be  TVD,  Q  must  be  sufficiently  large  (see  Eq.  (5)).  For  the  matrix 
viscosity  we  chose  Q  =  \A\  (modified  near  zero  eigenvalues)  while  for  the  scalar  viscosity  we 
chose  Q  =  u[A)I. 

For  the  CUSP  scheme  we  instead  choose  d  as  a  linear  combination  of  w  and  /.  We  shall 
only  consider  the  choice  for  the  state  vector  given  by 


Wh  =  (p  pu  pH) 

T 

(  ^  ] 

<  0  ^ 

f  =  u\  pu 

p 

[pH  ) 

This  choice  is  denoted  HCUSP  by  Jameson  [7].  The  first-order  accurate  CUSP  scheme  is 
defined  as 


dj+1/2  =  X^'C(U)J4.1  -  Wj]  +  -(fj+1  -  fj). 


(12) 


The  factor  c  is  included  so  that  i/  is  dimensionless.  We  thus  consider  only  scalar  parameters 
instead  of  a  matrix  coefBcient,  but  we  have  two  free  parameters,  v  and  The  scheme  is  total 
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enthalpy  preserving*  By  using  the  arithmetic  average,  u  =  +  Uj),  and  the  definition 


ac  =  i/c  + 


(13) 


one  can  rearrange  Eq*  (12)  to  obtain 

<^^+1/2  =  -Wj)  +  -  fpj)  +  -  Ujl 

Introducing  the  Roe  matrix  [17]  Arh  we  have  that  Jr  —  Jl  =  Arr(vjr  —  wr).  This  relation 
is  exact  if  Ari  is  computed  from  weighted  averages  of  the  left  and  the  right  states.  Then  the 
first-order  dissipation  is 

(ij+l/2  =  +  VCI){WR  -  Wl)-  (14) 

We  see  from  this  formula  that  d  is  a  linear  function  of  A.  Recall  that  m  is  a  quadratic  function 
of  d,  by  the  Cayley-Hamilton  Theorem.  Hence,  it  is  not  possible  to  bound  d  by  \Al  Thus, 
the  CUSP  scheme  cannot  be  TVD. 

Assume  that  the  subscript  L  denotes  the  interior  point  inside  the  shock  zone,  R  is  the 
state  downstream  of  the  shock,  and  the  state  LR  is  subsonic.  Jameson  [7]  shows  that  the 
downstream  point  with  the  state  R  is  in  equilibrium  if 

i/c 

fR-  fh  A  ~  ^ 

Substituting  the  Roe  matrix  for  the  difference  in  /  into  Eq.  (15)  we  get 

{^Ari  +  {'^’R  ~  '^r)  —  (1®) 

Hence,  wr  -  wl\s  sn  eigenvector  of  and  -(i/c)/(l  +  /?)  is  the  corresponding  eigenvalue. 
The  eigenvalues  of  Arl  are  known  to  be  A+,  and  u.  If  A  is  an  eigenvalue  of  A,  then  using 
this  formula  for  vc  in  Eq.  (14)  we  have 

<^^•+1/2  =  \  [-■^  +  P{Arl  -  A/)]  {WR  -  wl). 

In  order  to  have  a  positive  diffusion  when  i6  >  0,  we  require  that  A  be  negative  (i.e.,  —  (i'c)/(l  + 
/3)  =  A  ).  Thus, 

-(1  +  /?)A".  (17) 

For  u  <  0  we  obtain  similarly 

i/c  =  (1  -  /?)A+.  (18) 

So,  we  have  reduced  our  two  free  parameters  to  one  free  parameter  by  demanding  a  one  point 
shock  profile.  More  generally,  Jameson  shows  that  one  obtains  a  shodr  profile  with  one  interior 
point  if  the  following  two  conditions  hold: 
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1.  When  the  flow  is  supersonic  through  the  shock  then  one  obtains  a  totally  upwind  flux, 
n.  The  artificial  dissipation  Q  satisfies  a  generalized  eigenvalue  problem 


(^RL  “  OchaQra)  (^•’R  ~  0 


at  the  exit  fixjm  the  shock. 


The  second  condition  is  satisfied  by  both  the  matrix  viscosity  and  the  CUSP  scheme;  however, 
the  scalar  viscosity  does  not  satisfy  the  first  condition.  We  again  note  that  the  TVD  condition 
Q  >  |.A|  is  satisfied  by  the  scalar  and  matrix  viscosities  but  not  by  CUSP. 

What  remains  to  be  done  is  to  choose  suitable  functions  for  P  and  vc  which  satisfy  the 
above  requirements.  Jameson’s  choice  for  which  is  based  upon  the  eigenvalues  corresponding 
to  the  acoustic  waves,  is  given  by 


P  = 


-hmax 
—  max 


0, 

0, 


u  X~ 


U—X+ 


sgn(M) 


if  0  <  M  <  1 
if  -1  <  M  <  0 
if  \M\  >  1. 


(19) 


The  cutolfe,  /?  >  0  for  u  >  0  and  /?  <  0  for  u  <  0,  ensure  that  the  pressure  terms  are  discretized 
centrally  for  small  Mach  numbers.  Shock  capturing  with  one  interior  point  is  obtained  by  taking 


|u|  if  P=0 

-(1  +  /3)A-  if  /3>0,  0<M<1 
(1  -  P)X+  if  P<0,  -1  <  M  <  0 
0  if  |M|  >  1. 


The  dissipation  coeflEicients  are  to  be  computed  from  Roe-averaged  quantities.  They  provide 
full  upwinding  for  supersonic  flow  {P  —  sgn(M),  v  —  0).  The  choice  of  vc  —  |u|  for  /3  =  0  yields 
a  continuous  dissipation  coefficient  in  the  subsonic  region,  and  it  does  not  smear  slip  lines  with 
|iil  close  to  zero.  This  makes  the  CUSP  formulation  attractive  for  viscous  flow  calculations  with 
boundary  layers.  However,  viscous  flows  are  usually  discretized  by  using  cells  with  large  aspect 
ratio.  It  is  well  known  that  this  situation  requires  larger  dissipation  scaling  in  the  direction  of 
the  leeg  cell  sides  than  given  by  lu|.  We  redefine  the  dissipation  coefficients  in  the  individual 
coordinate  directions.  For  the  ^-direction  we  have 


vc^ 


max(|u|,5cr  ) 
-(l+i0)A- 


if 

if 


(l-i(?)A+  if 


if 


p  =  o 
P  >0  and 
0<M<1 
P  <0  and 
-1  <M<0 
\M\  >  1. 


(20) 


where  and  r  axe  functions  of  the  spectral  radii  in  the  ^  and  rj  directions  (A^  and  At^),  and 
they  are  defined  as  follows: 

r  =  =  niax(r,  1),  r  =  niin(r,  1). 


The  dissipation  coeSicient  in  the  T]!-direction  is  defined  correspondingly. 
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Simplified  Scheme 

Several  modifications  of  the  CUSP  scheme  have  been  in.  use  so  far.  Based  upon  the  Wh  = 
{p  pu  pH)^  system  the  dissipation  coefiicients  presented  in  Refs.  [7]  and  [23]  are  as  follows: 


a 


r  |M|  if 

U{.+^)  if 

+  raaxfo,^ 
sgn(M) 


\M\  >  e 
\M\  <  e, 

if  0  <  M  <  1 
if  -1  <  M  <  0 
if  \M\  >  1. 


(21) 


This  choice  does  not  allow  exact  shock  capturing  because  Eqs.  (17)  and  (18)  are  not  satisfied. 
Furthermore,  Roe  averaging  has  been  replaced  by  arithmetic  averaging  in  Ref,  [7]  and 
by  u  -  c,u-\-  c,  respectively  This  simplification  saves  a  few  square  roots  in  the  coding  of  the 
dissipative  flux.  This  then  becomes 


r  \M\  if  |Ml>e 
“  =  \  i(e  +  ^)  if  lM|<e, 

(22) 

[  max(0,2M-l)  if  0<M<1 
13  =  \  min(0,2M  +  l)  if  -1  <  M  <  0 
sgn(M)  if  \M\  >  1. 

Higher  Order  Scheme 

Having  determined  vc  and  /?,  we  see  from  Eq.  (12)  that  the  scheme  is  completely  defined  in 
terms  of  w  and  /,  Recalling  from  previous  sections,  we  now  need  to  combine  a  first-order 
artificial  dissipation  with  a  hi^-order  dissipation.  As  previously  the  high-order  dissipation  is 
used  in  the  smooth  re^ons  while  the  low-order  artificial  viscosity  is  used  in  the  neighborhood 
of  shocks.  To  allow  changing  from  one  type  to  the  other  we  can  either  use  the  JST  switch  of 
Eq.  (4)  or  the  SLIP  scheme.  Formula  (12),  as  given,  is  only  first-order  accurate,  as  it  depends 
only  on  dj+1/2  =  wj+i  -  Wj,  and  so  the  complete  artificial  viscosity  behaves  like  a  second 
difference.  In  order  to  improve  this  situation  we  need  to  use  a  MUSCL  or  SLIP  approach 
to  get  the  first  difference  to  higher  order  accuracy.  Since  the  fluxes  in  the  CUSP  scheme 
have  different  slopes  on  the  two  sides  of  the  sonic  line,  the  SLIP  construction,  which  relies 
on  successive  differences,  can  cause  a  loss  of  smoothness  at  the  sonic  line.  This  difficulty  can 
be  avoided  by  using  the  MUSCL  approach  to  obtain  the  states  corresponding  to  higher  order 
accuracy.  To  impose  monotonicity  we  can  apply  the  limiter  of  the  SLIP  scheme.  In  particular, 
we  can  replace  Eq.  (12)  by 

<fj+l/2  =  f  -  f{wL)]  , 

where  wr,wl  are  given  by  Eq.  (8).  This  procedure  was  followed  throughout  the  numerical 
examples  shown  below.  Application  of  Eq.  (8)  to  the  =  {p  pu  variables  still  allows 
total  enthalpy  to  be  preserved  in  the  higher  order  scheme. 
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Analysis  of  CUSP  Scheme 

The  eigenvalues  of  /SA/jl  +  vcl  in  Eq.  (14)  are  juic,  jU2C,  and  nzc.  Using  the  simplifications  of 
Eq.  (22)  the  eigenvalues  are: 

m  =  |M|, 

if  \M\  <  \ 
if  \<M<\ 
if  \M\  >  1, 

if  |M|  <  \ 

if  5  <  M  <  1 
if  \M\  >  1 

We  note  that  for  |M|  <  2  aU  three  eigenvalues  of  the  artificial  viscosity  are  equal,  and  so  we 
have  the  equivalent  of  a  scalar  viscosity.  The  scalar  viscosity  now  scales  with  \M\c  rather  than 
{\M\  +  l)c  for  the  JST  scalar  viscosity.  This  is  more  similar  to  the  case  of  preconditioning 
where  all  the  eigenvalues  are  appKsximately  \M\c  for  low  speed  flow.  Hence,  we  expect  that 
the  CUSP  dissipation  should  work  properly  for  very  low  Mach  numbers  provided  the  central 
flux  terms  are  augmented  by  a  suitable  preconditioning  matrix. 

In  the  subsonic  range  where  /?  =  0,  all  of  the  versions  of  the  CUSP  scheme  do  not  satisfy 
Eqs.  (17)  and  (18),  which  are  necessary  for  shock  capturing.  That  is,  the  cell-face  Mach 
numbers  in  the  shock  structure  have  to  be  larger  than  about  0.5  in  order  to  avoid  post-shock 
oscillations.  The  motivation  to  design  i/c  =  |ti|  for  /3  =  0  has  already  been  discussed  previously 
and  is  not  repeated  here.  However,  the  choice  of  the  function  for  as  given  in  Eq.  (19),  is 
not  necessarily  optimal.  For  example,  /3  =  max{0,{u  +  \\~)/{u  —  A“))  would  allow  shock 
capturing  for  Mach  numbers  down  to  about  1/3,  but  the  subsonic  dissipation  would  be  twice 
as  large,  vc  =  2lu|  for  /3  =  0.  Nevertheless,  our  own  experience  gained  from  a  number  of 
numerical  applications  suggests  that  there  is  no  need  for  further  modifications  of 

It  is  rather  difficult  to  compare  the  effect  of  the  parameter  of  the  JST  and  the  combined 

CUSP / SLIP  schemes  since  these  schemes  also  include  eigenvalue  information  which  is  not  the 
same.  To  isolate  the  effect  we  consider  a  bw  Mach  number  flow  with  preconditioning  (see 
Ref.  [25]  for  details).  Now  both  switches  are  based  on  the  convective  eigenvalue  u.  A  typical 
value  for  the  JST  scheme  is  However,  for  a  aspect  ratio  of  one  the  Martinelli  scaling 

[13]  adds  another  factor  of  two.  The  parameter  a  =  0  in  the  preconditioning  adds  an  additional 
factor  of  approximately  2.6.  Hence,  the  effective  constant  multiplying  the  fourth  difference  is 
about  which  is  somewhat  smaller  than  the  \  used  with  the  original  CUSP  scheme  with  the 
SLIP  formulation.  For  transonic  flows  it  is  more  difficult  to  compare  the  levels  of  dissipation. 
However,  it  seems  that  the  original  SLIP  scheme  yields  too  high  a  viscosity  level  and  so 
should  be  reduced  to  less  than  Numerical  computations  demonstrate  the  improved  accuracy 
(though  slower  convergence)  for  standard  transonic  turbulent  flows  when  is  reduced. 

Numerical  Results 

In  this  section  we  assess  the  accuracy  and  shodi  capturing  capabilities  of  the  HCUSP  scheme. 
Comparisons  are  made  between  the  HCUSP  and  MATD  schemes.  The  commonly  used  scalar 
dissipation  scheme  is  also  included  in  some  of  the  comparisons.  In  so  doing  one  can  clearly  see 


\M\ 

P2  =  <  OH-/? 

iM-bl 

'  \M\ 

HZ  =  I  a -13 
\M-1 
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HCUSP 
exact  dissipation 


0-5  x/C 


HCUSP  ^ 

approx,  dissipation 


0.5  x/c 


Influence  of  HCUSP  dissipation  coefficients  on  transonic  flow  over  NACA  0012  airfoU. 


Figure  2a:  Inviscid  pressure  distributions  computed  with  MATD  scheme  (NACA  0012  airfoil, 
M  =  0.80,  a  =  L25°), 

the  superiority  of  the  high-resolution  HCUSP  and  MATD  schemes  on  even  relatively  coarse 
meshes  (i,e,,  16  cells  in  the  boundary  layer  of  a  viscous  flow).  The  flow  problems  considered 
in  the  evaluation  of  these  numerical  diffusion  schemes  include  the  following:  1)  Inviscid  flow 
over  airfoils,  2)  laminar  flow  over  a  flat  plate,  3)  turbulent  flow  over  an  airfoil,  4)  inviscid 
and  viscous  hypersonic  flow  over  a  2-D  wedge.  The  computational  effort  and  convergence 
behavior  in  computing  these  solutions  is  given.  In  all  cases  a  five-stage  Runge-Kutta  scheme 
in  conjunction  with  the  convergence  acceleration  techniques  of  local  time  stepping,  implicit 
residual  smoothing,  and  multigrid  was  used. 

The  first  case  is  similar  to  the  application  published  in  Ref.  [7] .  Results  for  inviscid  transonic 
flow  over  the  NACA  0012  airfoil  are  shown  in  Fig.  1.  The  free-stream  Mach  number  is  0.8  and 
the  angle  of  attack  is  1.25°.  An  0-topology  mesh  of  160  x  32  cells  was  used  and  is  partially 
depicted  in  Fig.  1.  The  shock  waves  on  the  upper  and  lower  surfaces  are  captured  with  a  single 
interior  point.  The  differences  between  the  accurate  dissipation  coefficients  of  Eqs.  (19)  and 
(20)  on  the  one  hand  and  the  simplified  coefficients  of  Eq.  (22)  on  the  other  are  small.  Careful 
examination  reveals  that  the  exact  formulation  is  somewhat  more  accurate  on  the  upper  surface 
where  the  shock  is  stronger.  These  inviscid  solutions  were  obtained  with  a  4-level  multigrid 
and  sawtooth-type  cycle.  They  converged  at  a  rate  of  about  0.90  per  multigrid  cycle. 

In  Figs.  2  and  3  the  HCUSP  and  MATD  schemes  are  compared  for  this  case.  Solutions 
were  computed  on  three  successively  finer  C-topology  meshes.  The  coarsest  mesh  contained 
192  X  32  cells,  with  160  cells  on  the  airfoil,  and  for  each  sequential  mesh  the  number  of  cells 
in  each  coordinate  direction  was  doubled.  The  principal  differences  between  the  solutions 
occur  at  the  shock  waves.  Since  the  MATD  scheme  uses  a  pressure  switch  for  all  the  flow 
equations,  it  cannot  capture  a  shock  with  a  single  interior  point.  It  requires  three  interior 
points.  Nevertheless,  the  resolution  of  the  stronger  upper  surface  shock  is  nearly  the  same 
for  both  the  MATD  and  HCUSP  schemes  on  the  384  x  64  and  768  x  128  meshes.  With  the 
MATD  formulation  there  is  some  smearing  on  the  192  x  32  mesh.  As  clearly  evident  in  Fig.  3 
the  HCUSP  scheme  allows  a  sharp  definition  of  the  weak  lower  surface  shock  and  the  Zierep 
singularity  that  immediately  follows.  The  aerodynamic  coefficients  calculated  with  the  MATD 
scheme  on  the  384  x  64  mesh  essentially  agree  with  those  for  the  finest  grid.  The  lift  and 
drag  coefficients  determined  with  the  HCUSP  scheme  on  the  corresponding  meshes  are  not  the 
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Figure  2b:  Inviscid  pressure  distributions  computed  with  HCUSP  scheme  (NACA  0012  airfoil, 
M  =  0.80,  a  =  1.25°). 


same,  with  the  finest  grid  values  approaching  those  obtained  with  the  MATD  scheme.  These 
results  probably  reflect  the  higher  level  of  background  dissipation  with  the  HCUSP  scheme. 

As  an  initial  evaluation  of  the  dissipation  schemes  for  viscous  flows  we  consider  low-speed 
(Moo  =  0.15)  flow  over  a  flat  plate  at  zero  incidence.  For  this  flow  the  Reynolds  number  per 
unit  length  is  10®,  The  computational  domain  is  a  rectangle.  With  respect  to  the  leading  edge 
of  the  plate,  the  domain  extends  two  plate  lengths  upstream  and  one  plate  length  downstream. 
The  upper  boundary  is  four  plate  lengths  above  the  plate.  Solutions  were  computed  on  the 
same  domain  and  grids  used  in  Ref.  [23].  Starting  with  the  finest  mesh,  coarser  meshes  were 
determined  by  successively  eliminating  every  other  mesh  line.  The  finest  grid  consists  of 
512  X  128  cells,  with  384  cells  on  the  plate.  In  the  direction  y  normal  to  the  plate  the  grid  is 
spaced  uniformly  in  the  boundary-layer  coordinate  ij  (rj  =  y/ReJ'^),  where  x  is  the  coordinate 
parallel  to  the  surface,  and  Rex  is  the  Reynolds  number  based  on  distance  from  the  leading 
edge  of  the  plate).  Thus,  there  is  constant  resolution  of  the  boundary  layer  at  each  location 
along  the  plate.  Outside  the  boundary  layer  the  grid  is  stretched  exponentially.  In  order  to 
resolve  the  region  in  the  vicinity  of  the  stagnation  point,  the  grid  is  clustered  at  the  leading 
edge  of  the  plate.  At  the  surface  of  the  plate  no-slip  and  adiabatic  boundary  conditions  are 
enforced.  Along  the  boundary  upstream  of  the  leading  edge,  a  symmetry  condition  is  applied. 
Characteristic  type  boundary  conditions  are  used  at  the  upstream,  downstream,  and  upper 
boundaries. 

A  comparison  of  the  velocity  profile  at  X/L  =  0.82  computed  with  the  the  scalar,  matrix, 
and  HCUSP  dissipation  forms  is  displayed  in  Fig.  4.  Even  with  just  eight  points  in  the  boundary 
layer  (64  x  16  grid)  the  MATD  and  HCUSP  schemes  nearly  replicate  the  Blasius  solution.  As 
demonstrated  in  Ref.  [1]  scalar  dissipation  can  produce  serious  contamination.  With  the  scalar 
dissipation,  more  than  32  points  are  required  in  the  boundary  layer  to  obtain  a  grid  converged 
soiution.  For  the  MATD  and  HCUSP  schemes  the  variation  of  the  errors  (relative  to  the  Blasius 
solution)  in  the  calculated  skin  friction,  displacement  thickness,  and  momentum  thickness  are 
shown  in  Figs.  5a  and  5b.  The  standard  definitions  [18]  are  used  for  these  boundary-layer 
quantities.  The  errors  in  all  the  boundary-layer  parameters  are  quite  similar  for  the  high- 
resolution  schemes.  This  is  not  surprising  since  both  schemes  have  a  scaling  factor  that  vanishes 
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Figure  3a;  MATD  solutions  near  lower  surface  shock  (NACA  0012  airfoil,  M 
1.25°). 
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Figure  3b:  HCUSP  solutions  near  lower  surface  shock  (NACA  0012  airfoil,  M 
1.25 


=  Y*SQRT(Re,)/X  ^  ti  =  Y*SQRT(ReJ/X  _  t|  =  Y*SQRT(Re,)/X 


4a:  Tangential  and  transverse  velocity  profiles,  X/L  =  0.82,  64  x  16  grid. 


4b:  Tangential  and  transverse  velocity  profiles,  X/L  =  0.82,  128  x  32  grid. 


4c:  Tangential  and  transverse  velocity  profiles,  X/L  =  0.82,  256  x  64  grid. 


Figure  4:  Boundary-layer  profiles  on  flat  plate  with  M  =  0.15  and  Re  =  10^, 
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Figure  5a:  Comp  lie  results  Mth  Blasius  solution  {M  —  0.15  and  Re 

10^). 
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1 

Figur^ 5b:  Comparison  of  HCUSP  scheme  results  with  Blasius  solution  (M  =  0.15  and  Re 
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Figure  6a:  Comparison  of  pressure  distributions  with  SCALAR,  MATD,  and  HCUSP  schemes 
on  160  X  32  grid  (RAE  2822  airfoU,  M  =  0.73,  a  =  2.79°,  Re  =  6.5  x  10®). 


Figure  6b:  Comparison  of  skin-friction  distributions  with  SCALAR,  MATD,  and  HCUSP 
schemes  on  160  x  32  grid  (RAE  2822  mrfoil,  M  =  0.73,  a  =  2.79°,  Re  =  6.5  x  10®). 
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as  the  surface  is  approached. 

We  next  consider  transonic  flow  over  the  RAE  2822  airfoiL  The  free-stream  Mach  number 
is  0.73,  the  angle  of  attack  is  2.79°,  and  the  Reynolds  number  baaed  on  the  airfoil  chord  is 

6.5  X  10^.  Transition  from  laminar  to  turbulent  flow  is  fixed  at  the  3%  chord  location.  The 
C-type  grids  used  are  as  follows:  (1)  160  x  32  with  128  cells  on  the  airfoil,  (2)  320  x  64  with 
256  cells  on  the  airfoil,  and  (3)  640  x  128  with  512  cells  on  the  airfoil.  The  outer  boundary  is 
located  20  chords  from  the  airfoil.  The  normal  spacing  at  the  surface  of  the  640  x  128  mesh  is 

7.5  X  10  ^  chords.  At  the  leading  and  trailing  edges  of  the  airfoil  the  mesh  is  clustered,  giving 
tangential  spacings  of  1,17  x  10“^  and  1.86  x  10”^  chords,  respectively.  These  spacings  are 
roughly  doubled  with  each  mesh  coarsening.  To  determine  the  effect  of  further  mesh  refinement 
a  calculation  was  performed  with  the  MATD  scheme  on  a  1280  x  256  grid. 

In  Figs.  6-8  the  pressure  {Cp)  and  surface  skin-friction  (C7/)  distributions  computed  with 
the  different  dissipation  schemes  for  the  range  of  meshes  described  are  shown,  along  with  the 
experimental  data  [2].  As  in  the  inviscid  cases  the  primary  differences  in  the  solutions  occur 
at  the  shock  wave.  On  the  coarsest  grid  (160  x  32)  both  the  SCALAR  and  HCUSP  schemes 
produce  a  solution  with  the  shock  too  far  upstream.  This  is  an  unexpected  result  for  the 
HCUSP  scheme.  The  acceleration  of  the  flow  upstream  of  the  shock  is  underpredicted  relative 
to  the  finest  grid.  In  Ref.  [21]  the  adverse  effect  of  a  smooth  limiter  on  the  accuracy  of  the 
solution  in  the  vicinity  of  flow  transition,  and  thus  on  the  acceleration  of  the  flow  upstream 
of  the  shock,  is  demonstrated.  Therefore,  such  a  result  with  the  HCUSP  scheme  could  be 
a  consequence  of  the  smooth  limiter  being  used.  With  both  the  SCALAR  and  the  MATD 
schemes  a  nonphysical  increase  in  the  skin-friction  solution  on  the  upper  surface  appears  at 
the  trailing  edge  of  the  airfoil.  This  behavior  does  not  occur  in  the  solution  obtained  with  the 
HCUSP  scheme.  The  computed  aerodynamic  coefficients,  including  the  pressure  and  friction 
contributions  to  the  total  drag,  are  given  in  Table  1.  On  each  mesh  the  lift  and  drag  coefficients 
corresponding  to  the  solution  obtained  with  the  MATD  scheme  exhibit  the  closest  agreement 
with  the  1280  x  256  grid  values.  There  are  only  small  discrepancies  in  the  coefficients  associated 
with  the  MATD  and  HCUSP  schemes  on  the  640  x  128  grid. 

Convergence  behavior  for  the  HCUSP  and  MATD  schemes  is  displayed  in  Fig.  9.  Five  levels 
of  multigrid  were  used  for  this  case  and  either  50  or  70  cycles  were  executed  on  two  coarser 
meshes  in  order  to  obtain  an  initial  solution.  On  the  320  x  64  grid  the  average  rate  of  reduction 
of  the  residual  for  both  schemes  is  about  0.92.  Figure  9  also  shows  the  effect  of  the  modification 
given  by  Eq.  (20)  to  i/c  in  the  HCUSP  scheme.  The  convergence  is  improved  by  using  the  2-D 
formulation  for  the  dissipation  coefficient  i/c.  Note  that  convergence  with  C  =  0  was  possible 
for  this  transonic  case  but  not  for  the  hypersonic  case  presented  below. 

The  fourth  case  is  hypersonic  2-D  flow  over  a  blunt  wedge.  Figure  10  displays  the  solutions 
obtained  for  viscous  and  inviscid  flow  using  identical  meshes  of  64  x48  cells.  Physical  diffusion  is 
so  large  that  the  shock  profile  is  significantly  smeared  in  the  viscous  result.  For  inviscid  flow,  on 
the  other  hand,  we  obtain  perfect  capturing  with  a  single  interior  point  in  the  shock  structure  by 
using  the  formulation  of  Eqs.  (19)  and  (20).  Detailed  comparisons  of  the  hypersonic  wedge  flow 
solutions  yielded  by  the  CUSP  scheme  and  AUSM  have  been  presented  in  Ref.  [14].  The  shock 
capturing  capabilities  of  both  schemes  are  essentially  equal.  A  comparison  of  shock  profiles 
for  the  exact  and  the  simplified  coefficients  is  given  in  Fig.  11.  We  have  chosen  the  first-order 
scheme  in  order  to  address  the  pure  shock  capturing  capability  of  the  CUSP  scheme  without 
interference  from  the  limiter.  The  simplified  dissipation  coefficients  of  Eq.  (22)  produce  strong 
oscillations  at  the  shock,  even  though  there  is  substantial  physical  diffusion  present.  Hence  it  is 
concluded  that  an  accurate  implementation  of  dissipation  coefficients  is  a  must  for  hypersonic 
flows  with  strong  shocks. 
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Comparisons  of  computation  times  indicate  that  the  CUSP  scheme  needs  about  40%  more 
computer  time  than  the  basic  scalar  dissipation.  The  MATD  scheme  only  requires  about 
15%  additional  time.  MATD  requires  less  CPU  time  primarOy  because  it  needs  only  a  single 
evaluation  of  the  limiter  function.  Due  to  lower  inherent  dissipation,  computations  with  the 
CUSP  formulation  converge  somewhat  slower  for  transonic  flows  than  those  with  simple  scalar 
dissipation.  The  major  advantage  of  the  CUSP  approach  is  that  it  is  more  accurate  and 
more  robust  than  scalar  ■viscosity.  Our  numerical  tests  indicate  that  the  accuracy  of  the 
CUSP  scheme  is  intermediate  between  scalar  and  matrix  dissipation  for  transonic  flows.  For 
hjrpersonic  flows  it  seems  to  be  more  robust  than  the  matrix  viscosity  even  though  it  is  not 
TVD. 

Since  the  CUSP  scheme  is  implemented  through  artificial  dissipative  terms,  it  does  not 
have  to  be  applied  at  each  stage  of  the  Runge-Kutta  method.  In  particular,  the  diffusive  fluxes 
can  be  evaluated  only  at  the  first,  third,  and  fifth  sta^s  of  a  five-stage  method,  just  as  is 
typically  done  for  the  scalar  numerical  dissipation. 

Concluding  Remarks 

The  CUSP  scheme  has  been  studied  and  analyzed.  A  detail  comparison  has  been  made  between 
the  CUSP,  MATD,  and  scalar  dissipation  schemes.  For  transonic  inviscid  flows  the  CUSP 
scheme  allows  better  resolution  of  shock  waves,  since  they  are  captured  with  one  interior 
point.  However,  the  aerodynamic  quantities  such  as  lift  and  drag  obtained  with  the  CUSP 
scheme  are  not  as  accurate  on  coarser  meshes  (i.e.,  320  x  64  ceUs  or  less)  as  those  calculated 
with  the  MATD  scheme.  Both  the  CUSP  and  MATD  formulations  give  high  accuracy  in 
the  computation  of  high  Re  number  flat-plate  flow.  For  transonic  viscous  flows  and  coarser 
meshes  the  accuracy  in  aerodynamic  coeflicients  is  better  with  the  MATD  scheme  than  with 
the  CUSP  scheme.  This  loss  in  accuracy  with  the  CUSP  scheme  on  coarser  grids  appears  to  be 
a  consequence  of  limiter  activation  (i.e.,  reduction  to  first  order).  Convergence  of  the  HCUSP 
scheme  has  been  improved  by  introducing  the  aspect-ratio  scaling  factor  of  Eq.  (20). 

With  our  present  choice  of  HCUSP  dissipation  coefficients  it  has  been  shown  that  the 
resolution  of  strong  shock  waves  occurring  in  hypersonic  flows  is  possible  whereas  the  simplified 
coefficients  that  were  published  previously  failed. 

The  CUSP  scheme  requires  roughly  40%  more  computer  time  than  the  scalar  scheme,  while 
the  MATD  scheme  needs  about  15%  more  time.  Using  only  pressure  in  the  switches  for  MATD 
rather  than  using  a  different  measure  for  each  equation  reduces  the  robustness  of  the  algorithm 
but  requires  less  computer  time.  (Convergence  behavior  with  the  CUSP  and  MATD  schemes 
is  similar.  For  hypersonic  flows  the  CUSP  scheme  seems  to  be  more  robust  than  the  MATD 
scheme. 
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Figure  7a:  Comparison  of  pressure  distributions  with  SCALAR,  MATD,  and  HCUSP  schemes 
on  320  X  64  grid  (RAE  2822  airfoU,  M  =  0.73,  a  =  2.79°,  Re  =  6.5  x  10^), 


Figure  7b:  Comparison  of  skin-friction  distributions  with  SCALAR,  MATD,  and  HCUSP 
schemes  on  320  x  64  grid  (RAE  2822  sdrfoil,  M  =  0.73,  a  =  2.79°,  Re  =  6.5  x  10^). 
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Figure  8a:  Comparison  of  pressure  distributions  with  SCALAR,  MATD,  and  HCUSP  schemes 
on  640  X  128  grid  (RAE  2822  airfoU,  M  -  0.73,  a  =  2.79°,  Re  =  6.5  x  10^), 


Figure  8b:  Comparison  of  skin-friction  distributions  with  SCALAR,  MATD,  and  HCUSP 
schemes  on  64D  x  128  grid  (RAE  2822  airfoil,  M  =  0.73,  a  =  2.79°,  Re  =  6.5  x  10^). 
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Dissipation 

Scheme 

Grid 

Cl 

Cd 

Cdp 

Cof 

SCALAR 

160  X  32 

0.8172 

0.01728 

0.01275 

0.004532 

320  X  64 

0.8331 

0.01743 

0.01194 

0.005487 

640  X  128 

0.8532 

0.01782 

0.01225 

0.005574 

MATD 

160  X  32 

0.8304 

0.01818 

0.01251 

0.005662 

320  X  64 

0.8538 

0.01808 

0.01250 

0.005571 

640  X  128 

0.8597 

0.01799 

0.01246 

0.005535 

1280  X  256 

0.8611 

0.01800 

0.01246 

0.005544 

HCUSP 

160  X  32 

0.7987 

0.01926 

0.01367 

0.005594 

320  X  64 

0.8493 

0.01831 

0.01263 

0.005679 

640  X  128 

0.8592 

0.01803 

0.01245 

0.005585 

Table  1:  Lift  and  drag  coeflicients  for  turbulent  flow  over  RAE  2822  airfoil  (RAE  2822,  Moo  =  0.73, 
Q  =  2.79°,  i?ec  =  6,5x10®). 
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9b:  Effect  of  modified  i/c  on  HGUSP  scheme* 


Figure  9:  Convergence  histories  for  turbulent  flow  computations  (RAE  2822  airfoil,  Moo  =  0.73, 
a  =  2.79°,  Rec  =  6.5  x  10^), 


22 


blunted  wedge 
^  M=10,  a=0° 

\TyT^=5,  Re,=10000 

grid  64x48 


AM=0.5 


If 

llllllllllll.’l/ 

liSl^ 


grid  64x48 

L 

center  line  flow 


/  M  =9,  cx=0 
T/r^=5:Re,=10000 

first  order  scheme 


^o-o-o-ckJ^' 


ach  no. 


HCUSP 


jressure 


-135  -130  -125  -1 

X 


0.5 

10. 

0.0 

8 

M 

-0.5 

6 

-1.0 

4 

-1.5 

2 

[-2.0 

Q 

Mach  no. 


-  O.C 

HCUSP  : 
appr.  dissip.  -  -Q. 


/WiS 


ressure 


Figure  11:  Influence  of  HCUSP  dissipation  coefficients  on  hypersonic  flow  over  2-D  wedge. 
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